Trabajo Final
Métodos Estadísticos Avanzados en Ciencia de los Datos
Métodos Estadísticos Avanzados en Ciencia de los Datos
Pablo Saldarriaga-Aristizabal, Nicolás Prieto-Escobar, Obed Ríos-Ruiz
Semestre 2020-1
Repositorio
Introduccion
El modelamiento del comportamiento de diferentes variables es un tema que ha sido estudiado en sectores energéticos, industriales, económicos y financieros. De allí se comienza a apreciar tanto la importancia que tienen los datos hoy en día al igual que las técnicas utilizadas para su modelamiento. La estadística es una disciplina que se preocupa por la recolección, organización, interpretación y análisis de datos, que, según su aplicación puede traer un gran impacto en la industria al momento de la toma de decisiones. En particular, diferentes técnicas estadísticas han sido utilizadas en el sector financiero, las cuales permiten modelar comportamiento de los clientes, acciones, entre otras variables.
Diferentes industrias dentro de su funcionamiento, deben presentar ante la superintendencia información relacionando los gastos y ventas que presentaron anualmente. Además, se presume que del mercado colombiano, como en los procesos de las industria está esa interacción con todo el mercado, es posible pensar que exista una relación entre las diferentes variables macroeconomicas (ej. PIB, TRM, Balance Fiscal, Indice de Desempleo, etc.) y estos montos de costos y gastos de las empresas. Debido a la cantidad de información con la que se cuenta (información de costos y gastos para diferentes empresas en colombia entre los años 2016 y 2018), se sabe que no se cuenta con información suficiente para la construcción de un modelo por empresa que permita ver la relación existente entre las variables macroeconómicas y las variables asociadas a costos y gastos de ventas. Por lo anterior, es posible considerar un conjunto de datos como la consolidación de la información de todas las empresas junto con la información macroeconomica para los años en estudio, así buscando construir un modelo general para realizar la modelación de estas variables reportadas ante la superintendencia para un conjunto de empresas cuya industria sea similar.
Por lo tanto, para el conjunto de datos meniconado anteriormente, se buscará modelar la información de costos y gastos de venta a partir de las variables macroeconomicas disponibles, al igual que analizar si al realizar alguna transformación a dichas variables resulta relevante al momento de la creación del modelo. Se utilizarán modelos lineales, comenzando con la evaluación del modelo lineal general, hasta la aplicación de modelos de efectos mixtos. Esta última estructura de modelos, es bastante usado al momento de tener individuos que comparten la misma información pero tienen una salida diferente (en nuestro caso, todas las empresas comparten la misma información de las variables macroeconómicas), por lo que utilizar este tipo de modelos resulta de gran interés ya que permite modelar tanto efectos a nivel de individuo como agregando un efecto aleatorio.
Clasificación Industrial Seleccionada
El sector de la construcción es uno de los más relevantes en la economía colombiana. En nuestro contexto nacional el sector es considerado como uno de los más vitales para el desarrollo del país y representa uno de los más importantes rubros en materia de produccción interna componiendo cada año de 6 a 7 por ciento del producto interno bruto total y hasta un 7.1% del total de ocupados a nivel nacional. Dicho sector es caracterizado por sus fluctuaciones estacionarias fuertemente influenciadas por los planes de infraestructura de gran escala y los planes de gobierno. Respecto al último trimestre de 2019 tuvo un aumento de 3.4%, uno de los más altos a nivel de america latina.
Si bien las estimaciones globales para el año 2020 en Colombia para esta industria eran positivas, la coyuntura del COVID-19 ha de perturbar fuertemente el sector, afectando con alto impacto a los importadores de materiales y a la demanada frente a los retrasos en la ejecución de obras. Desde enero el sector de la producción de concreto ya estaba presentando caídas significativas de hasta un 8.3% respecto al año pasado en el mismo mes, por lo que se esperan peores resultados al cierre del segundo trimestre del año presente. Muchas compañías planearon incrementos en sus precios, con la esperanza de generar mayores ingresos, pero dichos planes han de ser postergados bajo la actual coyuntura. La demanda, el driver más relevante en la industria, claramente se ve desplazado efecto del impedimento de comercialización y la paralisis en el país bajo las medidas de cuarentena nacional. Un punto importante es que algunas compañías del sector podrán seguir con sus actividades producto de la inclusión de actividades de infraestructura como vital durante la crisis del COVID-19.
Frente a la incertidumbre que generan estos escenarios y como la dinámica particular de cada compañía que aporta al crecimiento de la industria total evoluciona en el tiempo se torna relevante construir análisis y modelos estadísticos robustos que respondan a las perturbaciones en el mercado y permitan obtener información valiosa para la toma de decisiones.
Consolidación de los datos (variables de salida y explicativas)
NOTA: toda esta parte inicial de preparación y consolidación de los datos se hizo en Python. Los programas para realizarla se encuentran en el repositorio. Este script de R toma solo el dataset final que fue resultado de toda la consolidación y preprocesamiento.
Los datos que de la variable de respuesta (costos y gastos de ventas de algunas empresas para los años 2016 a 2018) se extrajeron de: http://pie.supersociedades.gov.co.
Se conservaron solo las empresas cuya clasificación industrial uniforme internacional tenga la palabra “construcción” en ella y cuyos NIT estuvieran incluidos en el archivo de empresas de 2018. Luego de aplicar este filtro, se eliminaron también todas aquellas empresas que en algún año entre 2015 y 2018 tuvieran algún valor vacío (NaN) ya sea en sus costos de venta o en sus gastos de venta, dejando en total 33 empresas distintas en el dataset. A pesar de que solo se analizan en este reporte los años 2016 a 2018, fue necesario también descargar 2015 para poder calcular la variable en diferencias para 2016. Para calcular las variables de costos y gastos de ventas en diferencia porcentual, se usó su variación porcentual entre un año y otro, la cual se calculo con esta fórmula:
Costo_de_ventas_dif(t) = (Costo_de_ventas(t) – Costo_de_ventas(t-1))/Costo_de_ventas(t)
Donde t era el año. Así, se obtuvieron Costos_de_ventas_dif para los años 2016, 2017 y 2018.
NOTA: A pesar de que en este trabajo los modelos que usamos son solo para el costo de ventas, aplicamos la misma transformación para el gasto de ventas para los análisis descriptivos.
Los valores para las variables macroeconómicas (variables explicativas) se extrajeron del Banco de la República y de la DANE. Se descargaron archivos de Excel para cada una de las variables macro. Algunas de estas variables macro ya venían como variación porcentual entre un año y el otro. Sin embargo, otras no venían de esta manera. A continuación, detallamos la forma final en la que quedaría cada una de las variables (luego de nuestro preprocesamiento interno en scripts de Python).
PIB: variación porcentual del PIB total entre un año y otro
TRM: variación porcentual de la TRM promedio entre un año y otro
Balance en Cuenta Corriente: variación porcentual del balance de cuenta corriente entre un año y otro
Balance Fiscal: variación porcentual del balance fiscal entre un año y otro
Tasa de Intervención: tasa de intervención promedio del año (como ya era porcentaje, se prefirió no usar variación)
Desempleo: tasa de desempleo promedio durante el año (como ya era porcentaje, se prefirió no usar variación)
Inflación: tasa de inflación anual
Nota: se menciona que, ya posteriormente en este reporte de R, se aplica una transformación de log(1+X) para todas las variables (tanto para las de salida como para las de entrada). Esto es para poder aportarle interpretabilidad a los modelos de regresión (los coeficientes de las regresiones representarán los cambios porcentuales en la variable de salida que se tendrán al cambiar en porcentaje las variables explicativas).
Lectura de los datos
all_data = read.csv("PreparacionDatos/Datos_completos.csv", encoding = "UTF-8")
NIT = all_data['NIT']
data = all_data[, -c(1,2)]Análisis descriptivo
Damos una visualizacion inicial al conjunto de datos con el que se va a trabajar:
| Costo.de.ventas | Gastos.de.ventas | Costo.de.ventas_dif | Gastos.de.ventas_dif | TRM | PIB | Desempleo | Inflacion | Tasa_Intervencion | Balance_CC | Balance_Fiscal |
|---|---|---|---|---|---|---|---|---|---|---|
| $60,626,083 | $13,978,834 | -0.10 | 0.14 | 11% | 7% | 9% | 6% | 7% | 35% | 27% |
| $3,785,517 | $60,608 | 0.07 | 0.33 | 11% | 7% | 9% | 6% | 7% | 35% | 27% |
| $39,634,059 | $5,645,885 | 0.34 | 0.69 | 11% | 7% | 9% | 6% | 7% | 35% | 27% |
| $24,304,534 | $8,977,594 | -0.05 | -0.04 | 11% | 7% | 9% | 6% | 7% | 35% | 27% |
| $31,280,967 | $6,702,215 | -0.05 | -0.18 | 11% | 7% | 9% | 6% | 7% | 35% | 27% |
| $11,886,427 | $1,050,755 | 0.13 | 0.03 | 11% | 7% | 9% | 6% | 7% | 35% | 27% |
En la siguiente tabla, vemos información descriptiva de las variables que se considerarán en la etapa de modelamiento del trabajo, vemos así que contamos con un total de 11 variables, cada una con magnitudes diferentes. Vemos que algunos indicadores tienen valores entre 0 y 1 tal como lo es el desempleo, mientras que otras variables representan dinero, tal como la trm, el monto de gastos y ventas. Así, tenemos una idea inicial de las características del conjunto de datos, por lo que se aplicarán lo métodos correspondientes cuando sea necesario, si lo que se realizará es sensible a la escala de los datos.
| Costo.de.ventas | Gastos.de.ventas | Costo.de.ventas_dif | Gastos.de.ventas_dif | TRM | PIB | Desempleo | Inflacion | Tasa_Intervencion | Balance_CC | Balance_Fiscal | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| nbr.val | 9.900000e+01 | 9.900000e+01 | 99.0000000 | 99.0000000 | 99.0000000 | 99.0000000 | 99.0000000 | 99.0000000 | 99.0000000 | 99.0000000 | 99.0000000 |
| nbr.null | 0.000000e+00 | 0.000000e+00 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 |
| nbr.na | 0.000000e+00 | 0.000000e+00 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 |
| min | 5.913660e+05 | 3.436800e+04 | -0.9289315 | -0.7898731 | -0.0334936 | 0.0656288 | 0.0922332 | 0.0318000 | 0.0435185 | -0.2740147 | -0.2287460 |
| max | 2.119141e+08 | 3.541196e+07 | 1.0883866 | 6.1607837 | 0.1117616 | 0.0734318 | 0.0968134 | 0.0575000 | 0.0710204 | 0.3516454 | 0.2726036 |
| range | 2.113227e+08 | 3.537760e+07 | 2.0173181 | 6.9506568 | 0.1452552 | 0.0078030 | 0.0045802 | 0.0257000 | 0.0275019 | 0.6256601 | 0.5013496 |
| sum | 4.755036e+09 | 6.487101e+08 | -3.7634585 | 6.2073859 | 2.6432290 | 6.9358217 | 9.3332892 | 4.2966000 | 5.8012057 | 7.4840537 | -0.5776688 |
| median | 3.527843e+07 | 5.204078e+06 | -0.0254994 | 0.0107538 | 0.0018298 | 0.0711158 | 0.0937804 | 0.0409000 | 0.0612552 | 0.1491588 | -0.0613627 |
| mean | 4.803067e+07 | 6.552627e+06 | -0.0380147 | 0.0627009 | 0.0266993 | 0.0700588 | 0.0942756 | 0.0434000 | 0.0585980 | 0.0755965 | -0.0058350 |
| SE.mean | 4.567598e+06 | 6.609484e+05 | 0.0279735 | 0.0665708 | 0.0062481 | 0.0003305 | 0.0001922 | 0.0010748 | 0.0011499 | 0.0263314 | 0.0210523 |
| CI.mean.0.95 | 9.064249e+06 | 1.311631e+06 | 0.0555124 | 0.1321077 | 0.0123991 | 0.0006559 | 0.0003814 | 0.0021329 | 0.0022820 | 0.0522538 | 0.0417776 |
| var | 2.065432e+15 | 4.324843e+13 | 0.0774690 | 0.4387361 | 0.0038648 | 0.0000108 | 0.0000037 | 0.0001144 | 0.0001309 | 0.0686408 | 0.0438768 |
| std.dev | 4.544702e+07 | 6.576354e+06 | 0.2783325 | 0.6623716 | 0.0621675 | 0.0032887 | 0.0019121 | 0.0106940 | 0.0114417 | 0.2619939 | 0.2094678 |
| coef.var | 9.462084e-01 | 1.003621e+00 | -7.3217010 | 10.5639935 | 2.3284333 | 0.0469421 | 0.0202817 | 0.2464057 | 0.1952567 | 3.4656883 | -35.8982759 |
Además de la información descriptiva presentados en la tabla anterior, podemos ver para cada una de las variables, su distribución de forma visual con la ayuda de la creación de histogramas de frecuencia.
Vamos a obtener la informacion relacionada a las medidas de centralidad y dispersion del conjunto de datos. Inicialmente, presentamos información del vector de medias y medianas que describen la centralidad del conjunto de datos. Posteriormente, consideramos las matrices de Covarianza y Correlación para tener una intuición de la variabilidad de la información que consideramos.
Vector de Medias y Medianas
|
|
Análisis desde variables absolutas
Boxplot
par(mfrow=c(1,2))
bxplot_costos = boxplot(Costo.de.ventas~Year, data = all_data)
bxplot_gastos = boxplot(Gastos.de.ventas~Year, data = all_data)
### Outliers para costos
nits_tab = t(all_data[(all_data$Costo.de.ventas %in% bxplot_costos$out),]$NIT)
rownames(nits_tab) <- c("NITs")
kable(nits_tab)%>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
column_spec(1,bold=T)| NITs | 860009694 | 900378893 | 860009694 | 900378893 | 860009694 | 900378893 |
### Outliers para gastos
nits_tab = t(all_data[(all_data$Gastos.de.ventas %in% bxplot_gastos$out),]$NIT)
rownames(nits_tab) <- c("NITs")
kable(nits_tab)%>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
column_spec(1,bold=T)| NITs | 801002644 | 801002644 | 800015615 | 801002644 |
Análisis desde variables diferenciales
Boxplot
par(mfrow=c(1,2))
bxplot_costos_dif = boxplot(Costo.de.ventas_dif~Year, data = all_data)
bxplot_gastos_dif = boxplot(Gastos.de.ventas_dif~Year, data = all_data)
### Outliers para costos dif
nits_tab = t(all_data[(all_data$Costo.de.ventas_dif %in% bxplot_costos_dif$out),]$NIT)
rownames(nits_tab) <- c("NITs")
kable(nits_tab)%>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
column_spec(1,bold=T)| NITs | 800236890 | 890909034 | 900437650 | 800236890 | 890909034 | 900389088 | 800236890 |
### Outliers para gastos dif
nits_tab = t(all_data[(all_data$Gastos.de.ventas_dif %in% bxplot_gastos_dif$out),]$NIT)
rownames(nits_tab) <- c("NITs")
kable(nits_tab)%>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
column_spec(1,bold=T)| NITs | 800081030 | 890904459 | 900389088 | 800045720 | 800236890 | 890311366 | 890909034 | 900234565 | 900364670 | 900378893 |
Análisis:
Al observar la distribución del comportamiento de las variables de ventas y gastos para los 3 años en consideración, podemos ver que no hay mucha variación entre el comportamiento de estas a lo largo del tiempo, vemos que en particular, tanto el monto de costo y gasto promedio (mediana) es constante para los 3 años, aunque se nota un aumento en los costos en el año 2018 para el caso de algunas empresas. Ademas, vemos para que caso de los Costos, las mismas 2 empresas (relacionadas a los NITs 860009694 y 900378893) fueron las que presentaron un maypr valor de costos en los 3 periodos, por lo que el boxplot los considera registros outliers. Mirando la distribución de los gatos, notamos que una empresa tuvo los gastos mas altos en los 3 periodos (relacionadas al NIT 801002644) pero adicional, en el último año, la empresa con el NIT 800015615 tuvo un aumento en sus gastos, por lo que es considerado también como un outlier por el boxplot.
Pasando a la distribución de las variables de diferencia relativa de costos y gastos, vemos para ambos casos los registros outliers son aquellos que tienen un valor negativo en su diferencia o un incremento porcentual muy alto. Podemos ver que para el caso del 2016 en la variable costos, las empresas outliers fueron las correspondiente a los NITs 800236890, 890909034, 900437650, de estas, la única que no presenta una diferencia negativa fue la primera (800236890), por lo que puede indicar un aumento considerable (raro) en el mercado en término de los costos de venta (practicamente el doble del año anterior), mientras que todas las otras empresas presentaron una disminución considerable en este mismo concepto. Para el caso del 2017, tres empresas fueron consideradas como outlier, entre ellas, dos que se consideraron en el periodo anterior (800236890, 890909034), las cuales ambas reportaron una reducción de aproximadamente 65% en sus costos respecto al año anterior, mientras que la otra empresa considerada como ergistro raro (900389088) tuvo un aumento en sus costos de un 88% (casi el doble de los costos del año anterior). Para el último periodo, solo se encuentra una de las empresas mencionadas en los dos periodos anteriores (800236890), la cual vuelve a reportar una disminución considerable en sus costos.
Ahora mirando la variable de diferencia de gastos, solamente 2 empresas tienen comportamientos alto, presentando un aumento de casi el doble de los gastos anuales, estas empresas corresponden a 800081030 y 890904459 para los años 2016 y 2017 respectivamente, en cambio para el 2018, se tuvo variabilidad en los gastos de ventas de 8 empresas que son considerados anormales tanto por sus aumentos como disminuciones, auqellas empresas son 900389088, 800045720, 800236890, 890311366, 890909034, 900234565, 900364670 y 900378893. En aquellas empresas que tuvieron aumentos considerables, reportan casi hasta un cuarto del aumento de gastos en comparación al año anterior, mientras que aquellas que reportan una disminución, vemos que sus gastos se disminuyeron hasta el doble del periodo pasado.
Asi, estos cambios en las industrias en sus costos y gastos de ventas, pueden estar directamente relacionados con algunos eventos que hayan hecho la gran variación en la materia prima requerida para sus productos (tipos de mezlcas de cemento por ejemplo), además que el sector seleccionado depende bastante de los planes de infraestructura que se tengan, por lo que también es factible que en los periodos analizados, hayan ocurrido cambios en estos planes dentro de la contratación de cada empresa, y sean estas las causas que nos lleven a ver diferencias tan grandes de periodo a periodo (tanto positivas como negativas)
Matriz de Covarianzas
cov_data = cov(data)
kable(formatC(cov_data,format = "e", digits = 2)) %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
scroll_box("100%", height = "480px")| Costo.de.ventas | Gastos.de.ventas | Costo.de.ventas_dif | Gastos.de.ventas_dif | TRM | PIB | Desempleo | Inflacion | Tasa_Intervencion | Balance_CC | Balance_Fiscal | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Costo.de.ventas | 2.07e+15 | 1.35e+14 | 2.34e+06 | -1.56e+06 | 8.13e+04 | 5.00e+03 | -2.64e+02 | 6.97e+03 | 1.91e+03 | 2.99e+04 | 2.89e+05 |
| Gastos.de.ventas | 1.35e+14 | 4.32e+13 | 3.58e+05 | 1.04e+05 | 3.18e+03 | 3.62e+02 | 1.48e+02 | -3.37e+02 | -8.58e+02 | -2.08e+04 | 1.34e+04 |
| Costo.de.ventas_dif | 2.34e+06 | 3.58e+05 | 7.75e-02 | -3.65e-02 | 1.14e-03 | 9.75e-05 | 2.24e-05 | -3.04e-06 | -1.27e-04 | -3.20e-03 | 4.39e-03 |
| Gastos.de.ventas_dif | -1.56e+06 | 1.04e+05 | -3.65e-02 | 4.39e-01 | 8.23e-04 | 2.07e-04 | 1.46e-04 | -5.03e-04 | -8.59e-04 | -2.03e-02 | 4.91e-03 |
| TRM | 8.13e+04 | 3.18e+03 | 1.14e-03 | 8.23e-04 | 3.86e-03 | 1.78e-04 | -6.97e-05 | 5.52e-04 | 4.28e-04 | 9.34e-03 | 1.30e-02 |
| PIB | 5.00e+03 | 3.62e+02 | 9.75e-05 | 2.07e-04 | 1.78e-04 | 1.08e-05 | -6.78e-07 | 1.56e-05 | 4.78e-06 | 7.94e-05 | 6.29e-04 |
| Desempleo | -2.64e+02 | 1.48e+02 | 2.24e-05 | 1.46e-04 | -6.97e-05 | -6.78e-07 | 3.66e-06 | -1.92e-05 | -2.19e-05 | -5.01e-04 | -2.02e-04 |
| Inflacion | 6.97e+03 | -3.37e+02 | -3.04e-06 | -5.03e-04 | 5.52e-04 | 1.56e-05 | -1.92e-05 | 1.14e-04 | 1.16e-04 | 2.61e-03 | 1.73e-03 |
| Tasa_Intervencion | 1.91e+03 | -8.58e+02 | -1.27e-04 | -8.59e-04 | 4.28e-04 | 4.78e-06 | -2.19e-05 | 1.16e-04 | 1.31e-04 | 3.00e-03 | 1.25e-03 |
| Balance_CC | 2.99e+04 | -2.08e+04 | -3.20e-03 | -2.03e-02 | 9.34e-03 | 7.94e-05 | -5.01e-04 | 2.61e-03 | 3.00e-03 | 6.86e-02 | 2.69e-02 |
| Balance_Fiscal | 2.89e+05 | 1.34e+04 | 4.39e-03 | 4.91e-03 | 1.30e-02 | 6.29e-04 | -2.02e-04 | 1.73e-03 | 1.25e-03 | 2.69e-02 | 4.39e-02 |
Matriz de Correlación
cor_data = cor(data)
kable(round(cor_data,2)) %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
scroll_box("100%", height = "480px")| Costo.de.ventas | Gastos.de.ventas | Costo.de.ventas_dif | Gastos.de.ventas_dif | TRM | PIB | Desempleo | Inflacion | Tasa_Intervencion | Balance_CC | Balance_Fiscal | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Costo.de.ventas | 1.00 | 0.45 | 0.18 | -0.05 | 0.03 | 0.03 | 0.00 | 0.01 | 0.00 | 0.00 | 0.03 |
| Gastos.de.ventas | 0.45 | 1.00 | 0.20 | 0.02 | 0.01 | 0.02 | 0.01 | 0.00 | -0.01 | -0.01 | 0.01 |
| Costo.de.ventas_dif | 0.18 | 0.20 | 1.00 | -0.20 | 0.07 | 0.11 | 0.04 | 0.00 | -0.04 | -0.04 | 0.08 |
| Gastos.de.ventas_dif | -0.05 | 0.02 | -0.20 | 1.00 | 0.02 | 0.10 | 0.12 | -0.07 | -0.11 | -0.12 | 0.04 |
| TRM | 0.03 | 0.01 | 0.07 | 0.02 | 1.00 | 0.87 | -0.59 | 0.83 | 0.60 | 0.57 | 1.00 |
| PIB | 0.03 | 0.02 | 0.11 | 0.10 | 0.87 | 1.00 | -0.11 | 0.44 | 0.13 | 0.09 | 0.91 |
| Desempleo | 0.00 | 0.01 | 0.04 | 0.12 | -0.59 | -0.11 | 1.00 | -0.94 | -1.00 | -1.00 | -0.50 |
| Inflacion | 0.01 | 0.00 | 0.00 | -0.07 | 0.83 | 0.44 | -0.94 | 1.00 | 0.95 | 0.93 | 0.77 |
| Tasa_Intervencion | 0.00 | -0.01 | -0.04 | -0.11 | 0.60 | 0.13 | -1.00 | 0.95 | 1.00 | 1.00 | 0.52 |
| Balance_CC | 0.00 | -0.01 | -0.04 | -0.12 | 0.57 | 0.09 | -1.00 | 0.93 | 1.00 | 1.00 | 0.49 |
| Balance_Fiscal | 0.03 | 0.01 | 0.08 | 0.04 | 1.00 | 0.91 | -0.50 | 0.77 | 0.52 | 0.49 | 1.00 |
corrplot(cor_data, method="circle")
Dado que en este trabajo se pretende realizar la modelación de las variables de costo y gasto de ventas (también considerando la variante en las diferencias de las variables por periodo), la matriz de correlación es un buen indicador para observar relaciones lineales existentes en las variables. En particular, vemos para el caso de la información de gastos, una “fuerte” correlación de las variables de costo, diferencia de costos y diferencia de ventas. Para este caso, es natural encontrar una correlacion positiva con dichas variables, pues es claro que costos y gastos están asociados entre sí. Igualmente, la variable de diferencia de gastos, fue calculada con la variable gastos, por lo tanto tiene sentido encontrar esta correlación. Por otro lado, tanto para la variable de ventas como gastos, no se ven correlaciones fuertes en relación a las variables macro-económicas, lo que nos indica que no existe una relación lineal entre esta variable y las demás. Por lo que una buena alternativa en este trabajo, será considerar transformaciones no lineales de las variables para encontrar una dependencia con la información de los costos de venta. Ahora, mirando las variables de diferencia de gastos y costos, vemos que existen dependencias lineales mayores con las variables macroeconómicas a diferencia de las variables originales, así vemos que la diferencia en costos tiene un grado de asociación con todas las variables macroeconomicas, mientras que los gastos presenta una relación lineal con el balance de cuenta corriente, por lo que se podría pensar que existen relaciones no lineales respecto a los otros indicadores económicos.
Visualizacion de la distribucion de los datos
pairs(data,diag.panel = panel.hist, lower.panel = panel.cor)
Visualizacion de los datos segun curvas de normalidad
### Funcion para encontrar los contornos
c_alpha = function(alpha, sigma, p){
res = (2*pi)^(-p/2)*(det(sigma))^(-1/2)*exp(-1/2*qchisq(1-alpha, df = p))
return(res)
}
grafica = function(data, name1, name2){
data_aux = data[c(name1,name2)]
names(data_aux) = c("y1","y2")
cov_data = cov(data_aux)
mean_data = mean=colMeans(data_aux)
min_value1 = min(data[name1])
max_value1 = max(data[name1])
min_value2 = min(data[name2])
max_value2 = max(data[name2])
n = 100
y1 = seq(min_value1, max_value1, length.out = n)
y2 = seq(min_value2, max_value2, length.out = n)
grid = expand.grid(y1,y2)
grid$Z<-apply(grid,1,dmvnorm,mean = mean_data,sigma=cov_data)
Z<-matrix(grid$Z,nrow=n,ncol=n)
contornos = sapply(c(0.01, 0.05, 0.1), c_alpha, sigma = cov_data, 2)
contour(y1,y2,Z,levels=contornos,labels=c("99%","95%","90%"),
las=1)
points(data_aux$y1,data_aux$y2)
grid()
title(main = "Contornos de distribucion normal", xlab = name1, ylab = name2)
}p1 = grafica(data, "Costo.de.ventas_dif", "PIB")
p2 = grafica(data, "Costo.de.ventas_dif", "TRM")
p3 = grafica(data, "Costo.de.ventas" , "Costo.de.ventas_dif")
p4 = grafica(data, "Gastos.de.ventas" , "Gastos.de.ventas_dif")
p4 = grafica(data, "Costo.de.ventas_dif" , "Gastos.de.ventas_dif")
Creacion de modelos
Para la variable costos de venta, adicional a la transformación que se realizó para obtener la diferencia relativa entre los periodos, en la etapa de modeloamiento, consideraremos para la variable salida la siguiente transformacion:
\[log(1+diff)\]
Adicionalmente, a los predictores se les realizó la misma transformación.
Agregar un comentario sobre la eleccion de las empresas para train y para test ### Pablo
data_aux = select(all_data,select = -starts_with('NIT'))
data_aux = select(data_aux,select = -starts_with('Year'))
### Transformacion de las variables respuesta
### Recordar para la interpretacion de los modelos:
### regresion log(y) = b*log(x) -> cambio en x, es un aumento en b% en y
### regresion log(y) = b*x -> cambio en x implica un aumento en 100*b puntos en y
### regresion y = b*log(x) -> cambio en x implica un cambio de (b/100)% en y
### regresion y = b*x -> cambio en x implica un aumento de b en y
### Transformacion 1: a nivel logaritmico
data_aux$Costo.de.ventas_dif = log(1 + all_data$Costo.de.ventas_dif)
data_aux$PIB = log(1 + all_data$PIB)
data_aux$TRM = log(1 + all_data$TRM)
data_aux$Desempleo = log(1 + all_data$Desempleo)
data_aux$Balance_CC = log(1 + all_data$Balance_CC)
data_aux$Balance_Fiscal = log(1 + all_data$Balance_Fiscal)
data_aux$Inflacion = log(1 + all_data$Inflacion)
data_aux$Tasa_Intervencion = log(1 + all_data$Tasa_Intervencion)empresas_train = c(800015615, 800045720,
800081030, 800112440,
800118660, 800157469,
800232356, 800236890,
801002644, 805012368,
806014553, 830030574,
830037495, 860009694,
860033653, 860050956,
890909034, 900173460,
900184722, 900204182,
900364670, 900378893)
empresas_test = c(830052054, 860030360,
860501682, 890117431,
890300012, 890311366,
890904459, 890929951,
900234565, 900389088,
900437650)train = all_data$NIT%in%empresas_train
test = all_data$NIT%in%empresas_test
NIT_train = NIT[train,]
NIT_test = NIT[test,]
data_train = data_aux[train,]
rownames(data_train) <- NULL
data_train_z = data_train
media_tr = attr(data_train_z, 'scaled:center')
stdev_tr = attr(data_train_z, 'scaled:scale')
data_train_z = as.data.frame(data_train_z)
data_test = data_aux[test,]
rownames(data_test) <- NULL
data_test_z = data_test
data_train_z$NIT = NIT_train
data_test_z$NIT = NIT_testDefinición de las funciones para calcular el R-squared y R-Squared Adjusted
r2_score <- function(x, y) summary(lm(y~x))$r.squared
adj_r2_score <- function(x, y) summary(lm(y~x))$adj.r.squared
mse_score <- function(x, y) mean((x - y)^2)Modelo para los Costos
Modelo Lineal General
step.model <- ols_step_both_p(lm('Costo.de.ventas_dif~PIB+TRM+Desempleo+Inflacion+Tasa_Intervencion+Balance_CC+Balance_Fiscal', data = data_train_z), pent = 0.1, prem = 0.1, details = TRUE)## Stepwise Selection Method
## ---------------------------
##
## Candidate Terms:
##
## 1. PIB
## 2. TRM
## 3. Desempleo
## 4. Inflacion
## 5. Tasa_Intervencion
## 6. Balance_CC
## 7. Balance_Fiscal
##
## We are selecting variables based on p value...
##
##
## Stepwise Selection: Step 1
##
## - PIB added
##
## Model Summary
## -----------------------------------------------------------------
## R 0.086 RMSE 0.457
## R-Squared 0.007 Coef. Var -429.959
## Adj. R-Squared -0.008 MSE 0.209
## Pred R-Squared -0.045 MAE 0.240
## -----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.100 1 0.100 0.478 0.4919
## Residual 13.367 64 0.209
## Total 13.467 65
## ------------------------------------------------------------------
##
## Parameter Estimates
## -----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -----------------------------------------------------------------------------------------
## (Intercept) -0.967 1.246 -0.776 0.441 -3.456 1.522
## PIB 12.709 18.384 0.086 0.691 0.492 -24.017 49.436
## -----------------------------------------------------------------------------------------## Warning in if (pvals[minp] <= pent) {: la condición tiene longitud > 1 y sólo el
## primer elemento será usado##
## No more variables to be added/removed.
##
##
## Final Model Output
## ------------------
##
## Model Summary
## -----------------------------------------------------------------
## R 0.086 RMSE 0.457
## R-Squared 0.007 Coef. Var -429.959
## Adj. R-Squared -0.008 MSE 0.209
## Pred R-Squared -0.045 MAE 0.240
## -----------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## ------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## ------------------------------------------------------------------
## Regression 0.100 1 0.100 0.478 0.4919
## Residual 13.367 64 0.209
## Total 13.467 65
## ------------------------------------------------------------------
##
## Parameter Estimates
## -----------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## -----------------------------------------------------------------------------------------
## (Intercept) -0.967 1.246 -0.776 0.441 -3.456 1.522
## PIB 12.709 18.384 0.086 0.691 0.492 -24.017 49.436
## -----------------------------------------------------------------------------------------summary(step.model)## Length Class Mode
## orders 1 -none- character
## method 1 -none- character
## steps 1 -none- numeric
## predictors 1 -none- character
## rsquare 1 -none- numeric
## aic 1 -none- numeric
## sbc 1 -none- numeric
## sbic 1 -none- numeric
## adjr 1 -none- numeric
## rmse 1 -none- numeric
## mallows_cp 1 -none- numeric
## indvar 7 -none- character
## betas 2 -none- numeric
## lbetas 1 -none- numeric
## pvalues 2 -none- numeric
## beta_pval 4 data.frame list
## model 12 lm listDe los resultados anteriores, vemos que el modelo lineal general con mejor desempeño es aquel que solo considera el PIB además de un intercepto. Pues vemos que en particular para el PIB, su error estandar se reduce considerablemente (un poco más de la mitad), además que vemos mejores valores al R-Squared Adjusted. Por lo tanto, procedemos a mirar el desempeño de un modelo entrenado considerando unicamente el PIB. ### ARREGLAR
mod_lin = lm('Costo.de.ventas_dif~PIB', data = data_train_z)En este caso el modelo está siendo exclusivamente dependiente de la dinámica el PIB, y se tiene que un cambio en dicha medida tiene un impacto positivo de hasta 12.70% sobre la variación en los costos de venta, lo cual es razonable, en particular porque al percibir mayor crecimiento en la economía global la industria, naturalmente, también se ve afectada al haber mayor producción y un mayor retorno de dinero; no obstante, las dinámicas de las compañías que componen la industria muchas veces se ven influenciada por el factor macro más global de manera diferentes, ello por su tamaño o tendencias, así que probablemente sea necesario modelar dichos efectos de una forma más robusta.
summary(mod_lin)##
## Call:
## lm(formula = "Costo.de.ventas_dif~PIB", data = data_train_z)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.55042 -0.01670 0.08328 0.19822 0.80263
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9668 1.2460 -0.776 0.441
## PIB 12.7095 18.3843 0.691 0.492
##
## Residual standard error: 0.457 on 64 degrees of freedom
## Multiple R-squared: 0.007412, Adjusted R-squared: -0.008097
## F-statistic: 0.4779 on 1 and 64 DF, p-value: 0.4919Desempeño del Modelo Entrenado
En el conjunto de entrenamiento:
preds = predict(mod_lin)
plot(lapply(preds, function(x) exp(x) - 1),
lapply(data_train_z$Costo.de.ventas_dif, function(x) exp(x) - 1),
xlab = 'Valor predicho', ylab = 'Valor real', main='Modelo Lineal General (a escala porcentual)')
abline(a=0, b=1, lwd = 2, col = 'red')
Como puede apreciarse en la gráfica anterior el comportamiento de las estimaciones es relativamente estática y no logra distribuir los puntos de información de forma dinámica. En particular, es claro como se obtuvo una segmentación en tres partes, lo cual es razonable y va en línea con la idea de crecimiento, estancamiento y decaida de los costos de venta de las compañías de la industria, es decir, si bien los valores predichos se encuentran desvíados de los reales, el sólo PIB y la propia variación de los costos permiten identificar tendencias en la industria.
La aplicación sobre el conjunto de prueba desplegado a continuación genera resultados muy similares a aquellos sobre el de entrenamiento. El ajuste es considerablemente bajo, lo cual es esperado, debido a la separación natural de las tendencias de las empresas y a la falta de la modelación del efecto del indicador sobre la medida de estudio.
preds_test = predict(mod_lin, newdata = data_test_z)
r2_model<-r2_score(preds_test, data_test_z$Costo.de.ventas_dif)
adj_r2_model<-adj_r2_score(preds_test, data_test_z$Costo.de.ventas_dif)
mse_model<-mse_score(preds_test, data_test_z$Costo.de.ventas_dif)
sub_tit = paste("R2", format(r2_model, digits=2, nsmall=2),
"; R2_adj", format(adj_r2_model, digits=2, nsmall=2),
"; MSE", format(mse_model, digits=2, nsmall=2),
sep = " ", collapse = NULL)
plot(lapply(preds_test, function(x) exp(x) - 1),
lapply(data_test_z$Costo.de.ventas_dif, function(x) exp(x) - 1),
xlab = 'Valor predicho', ylab = 'Valor Real', main='Modelo Lineal General (a escala porcentual)', sub=sub_tit)
abline(a=0, b=1, lwd = 2, col = 'red')
Graficos de Evaluacion de modelos
par(mfrow=c(3,2))
mod_plot = plot(mod_lin, which = c(1:6))
De los gráficos anteriores, en particular analizando el valor de la distancia de Cook para las observaciones, notamos que las observaciones 17, 30, 52 son candidatas a ser registros atípicos en el conjunto de datos, por lo que entrenaremos un nuevo modelo eliminando estos registros. Notar que los outliers corresponden a la empresa con NIT 890909034 para el 2016, y para la empresa 800236890 en los periodos del 2017 y 2018, información que se contrasta con los outliers obtenidos en el boxplot, lo cual están incluidos tanto mirando la variable de diferencia costos de venta.En vista de su comportamiento raro en magnitud según el boxplot, no es de extrañarse que estos registros hayan sido aquellos con la distancia de cook mas grande. Procederemos a entrenar un modelo removiendo estos registros
Modelos Lineales Generalizados (GLM)
Dentro de los modelos que fueron expuestos a lo largo del curso, se encuentra la familia de los modelos lineales generalizados, los cual nos permiten la creación de modelos teniendo en cuenta diferentes supuestos que se hacen sobre la variable respuesta que estamos considerando. En la elaboración de este trabajo, se está considerando realiza un modelo de predicción sobre transformaciones en la diferencia relativa de los costos de ventas que tienen las empresas, por lo que el dominio de la variable son los números reales (ya que puede tomar valores continuos tanto positivos como negativos, los cuales son considerados en las diferentes transformaciones que se plantean). Justo por la caracteristica de la variable respuesta que tenemos, no se considera pertinente realizar evaluacion de los modelos, esta conclusión está apoyada del análisis que se realiza de los diferentes modelos generalizados que podriamos considerar:
Poisson: Para este modelo, esperamos que la variable respuesta tenga la forma de conteos, no negativa. Por lo que la diferencia de los costos de venta no aplica para este modelo.
Logit: Este es un modelo logistico, el cual asume que la variable respuesta tendrá un comportamiento dentro del intervalo (0,1) el cual se utiliza para la predicción de la ocurrencia de un evento (ocurre o no ocurre), por lo que tampoco se ajusta para el modelamiento de la diferencia de costos de venta
Gamma y Gaussiana inversa: De las distribuciones de probabilidad de la función gamma y gaussiana inversa, se sabe que los valores de la variable pueden ser continuos, por lo que nos llevaría a pensar de forma inicial que alguno de estos se puede considerar, pero en vista que tienen otra testricción, y es que los valores son continuos positivos, la variable de diferencia de costo de ventas tampoco se podría modelar con este tipo de modelos, ya que puede tomar valores negativos
Gaussiana: Por la caracteristica de la distribución, sabemos que es posible modelar variables respuesta que su dominio sea los números reales, por lo que la variable de diferencia de costos de venta encaja en este tipo de modelo, ahora bien, esto es equivalente a considerar el modelo lineal general presentado anteriormente.
Por el análisis realizado, no se considera pertinente utilizar alguno de los otros modelos lineales generalizados para el modela miento de la variable transformada de diferencia relativa de costos de venta
Modelo de Efectos Mixtos
mod_me = lmer('Costo.de.ventas_dif~Desempleo+PIB+(0+Desempleo|NIT)+(PIB|NIT)', data = data_train_z)
summary(mod_me)## Linear mixed model fit by REML ['lmerMod']
## Formula: Costo.de.ventas_dif ~ Desempleo + PIB + (0 + Desempleo | NIT) +
## (PIB | NIT)
## Data: data_train_z
##
## REML criterion at convergence: 68.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2577 -0.0832 0.1859 0.3549 2.5374
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## NIT Desempleo 3.576169 1.89108
## NIT.1 (Intercept) 0.004193 0.06475
## PIB 0.566302 0.75253 -1.00
## Residual 0.182955 0.42773
## Number of obs: 66, groups: NIT, 22
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.8466 3.0990 -0.273
## Desempleo -1.2765 30.4740 -0.042
## PIB 12.6315 17.3072 0.730
##
## Correlation of Fixed Effects:
## (Intr) Desmpl
## Desempleo -0.926
## PIB -0.473 0.107Ahora, considerando el modelo de efectos mixtos, observando primero la componente de los efectos fijos, vemos que un incremento en un 1% del PIB anual del pais, refleja un aumento en un 12.63% en la diferencia de los costos de venta anuales de las empresas, mientras que al ver un cambio de un 1% en el desempleo, este disminuye la diferencia de los costos en un 1.2765%. Dichos resultados son completamente razonables en la realidad por el argumento previamente expuesto sobre el efecto del crecimiento del PIB en la variación de los costos al presentarse mayor producción que genera incrementos en los retornos de dinero, por otro lado, una variación negativa sobre el desempleo combinado con el PIB claramente afecta la industria al presentarse una mano de obra más barata lo que reduce los gastos y de forma indirecto los costos. Aquí es importante notar como la tasa de crecimiento potencial asociada al crecimiento de productividad debe, en teoría, mantenerse en línea con una reducción con la tasa de desempleo.
Ahora bien, considerando los efectos aleatorios, los cuales tendrán un efecto sobre la estimación en aquellas empresas que fueron consideradas en el entrenamiento (mirando al NIT), vemos que adicional a la componente fija, el cambio de 1% del PIB anual, que antes tenia un efecto en el aumento del 12.63% sobre la diferencia de los costos de venta, ahora se aumenta en un 0.56% adicional, es decir, para aquellas empresas que haya visto el modelo previamente, el aumento anual en los costos de venta corresponderá a un 13.198%. Ahora realizando el análisis de los efectos aleatorios para el desempleo, vemos que el desempleo tiene un incremento del 3.576% en la diferencia de los costos de venta anual por un aumento del 1% en el desempleo, por lo que para aquellas empresas consideradas al momento de entrenar el modelo, presentarán un aumento del 2.299669% en la diferencia de costos anuales por un cambio de un 1% en el desempleo, por lo que estas empresas pueden tener caracteristicas que evitan que el desempleo no tenga un imacto muy grande en sus costos, lo cual, una vez más, refleja la realidad de la industria de la construcción donde se presentan actores de magnitudes muy diferentes que pueden mitigar el efecto macro al poder mantener niveles de producción altos que les permita seguir generando imgresos sin tener que afectar a su capital humano.
preds = predict(mod_me)
plot(lapply(preds, function(x) exp(x) - 1),
lapply(data_train_z$Costo.de.ventas_dif, function(x) exp(x) - 1),
main='Modelo de Efectos Mixtos (a escala porcentual)', xlab = 'Valor predicho', ylab = 'Valor Real')
abline(a=0, b=1, lwd = 2, col = 'red')
preds = predict(mod_me, newdata = data_test_z,allow.new.levels=TRUE)
r2_model<-r2_score(preds, data_test_z$Costo.de.ventas_dif)
adj_r2_model<-adj_r2_score(preds, data_test_z$Costo.de.ventas_dif)
mse_model<-mse_score(preds, data_test_z$Costo.de.ventas_dif)
sub_tit = paste("R2", format(r2_model, digits=2, nsmall=2),
"; R2_adj", format(adj_r2_model, digits=2, nsmall=2),
"; MSE", format(mse_model, digits=2, nsmall=2),
sep = " ", collapse = NULL)
plot(lapply(preds, function(x) exp(x) - 1),
lapply(data_test_z$Costo.de.ventas_dif, function(x) exp(x) - 1),
main='Modelo de Efectos Mixtos (a escala porcentual)', xlab = 'Valor predicho', ylab = 'Valor Real', sub=sub_tit)
abline(a=0, b=1, lwd = 2, col = 'red')
Desafortunadamente en el caso de prueba no se logra distinguir la tendencia natural en el comportamiento de la diferencia de los costos de venta y la escala se encuentra muy distante entre los valores predichos y los reales, como se espera desde el gráfico anterior. En esta evaluación los valores predichos se ponderan para generar una tendencia casi constante y la razonabilidad una vez más se pone en duda, aunque no por completo, dado que éstos valores efectivamente pueden ser alcanzados con condiciones ligeramente diferentes.
Superficie poblacional del modelo
Miraremos la superficie poblacional del modelo, es decir, como cambia el valor predicho de una solución general desconcodia de la población (para la cual solo se le aplicarían los efectos fijos, y no los aleatorios) al variar los valores de las dos variables explicativas que inciden sobre ella (desempleo y PIB). Recordar que todas las variables en estos gráficos tienen ya la transformación del logaritmo antes explicada.
Desempleo = seq(min(data_train_z$Desempleo),max(data_train_z$Desempleo), length.out = 50)
PIB = seq(min(data_train_z$PIB),max(data_train_z$PIB), length.out = 50)
data_surface = expand.grid(Desempleo,PIB)
names(data_surface) = c('Desempleo','PIB')
data_surface$NIT = 1000
z = predict(mod_me, newdata = data_surface,allow.new.levels=TRUE)
z = matrix(z, nrow = 50, ncol = 50)
custom_txt <- paste0("Desempleo: ", data_surface$Desempleo,
"\nPIB: ", data_surface$PIB, # correct break syntax
"\nCosto_de_ventas_dif: ", z) %>%
matrix(50,50) # dim must match plotly's under-the-hood? matrix
### Toca rotar la figura para poder ver el plano
fig <- plot_ly(x = PIB, y = Desempleo, z = z, text = custom_txt, hoverinfo = "text") %>% add_surface(colorbar=list(title='Poblacional'))%>%layout(scene = list(xaxis = list(title = 'PIB'), yaxis = list(title = 'Desempleo'),zaxis = list(title = 'Costo_de_ventas_dif')))
figSuperficie del modelo para un individuo conocido
Para visualizar el funcionamiento del modelo de efectos mixtos, además de graficar la superficie poblacional general, graficaremos también la superficie de un par de invidiuos ya conocidos por el modelo. En este caso, como los individuos fueron usados por el modelo de efectos mixtos al entrenarse, se les aplica a dichos individuos un efecto aleatorio propio a cada uno de ellos (además del efecto fijo general que se aplica de la misma forma para todos). Por esta razón, veremos que las superficies de estos dos individuos son distintas entre sí y ademas son distintas de la superficie poblacional general.
Desempleo = seq(min(data_train_z$Desempleo),max(data_train_z$Desempleo), length.out = 50)
PIB = seq(min(data_train_z$PIB),max(data_train_z$PIB), length.out = 50)
data_surface = expand.grid(Desempleo,PIB)
names(data_surface) = c('Desempleo','PIB')
data_surface$NIT = empresas_train[1]
z2 = predict(mod_me, newdata = data_surface)
z2 = matrix(z2, nrow = 50, ncol = 50)
### Toca rotar la figura para poder ver el plano
data_surface$NIT = empresas_train[3]
z3 = predict(mod_me, newdata = data_surface)
z3 = matrix(z3, nrow = 50, ncol = 50)
fig %>% add_surface(z = ~z2, opacity = 0.98,colorscale = list(c(0,1),c("rgb(255,112,184)","rgb(120,0,64)")),colorbar=list(title=paste('NIT ', empresas_train[1]))) %>% add_surface(z = ~z3, opacity = 0.94,colorscale = list(c(0,1),c("rgb(0,0,0)","rgb(200,200,200)")), colorbar=list( title=paste('NIT ', empresas_train[1]))) %>% layout(title="Poblacional (azul y verde) vs dos Individuos (uno rosa y uno negro)")Conclusiones
Referencias
De donde sacamos los datos